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Abstract 

The Broad Histogram Method (BHM) allows one to determine the en- 
ergy degeneracy g(E), i.e. the energy spectrum of a given system, from the 
knowledge of the microcanonical averages < N up (E) > and < N dn (E) > of 
two macroscopic quantities N up and iV dn defined within the method. The 
fundamental BHM equation relating g(E) to the quoted averages is exact 
and completely general for any conceivable system. Thus, the only possible 
source of numerical inaccuracies resides on the measurement of the averages 
themselves. 

In this text, we introduce a Monte Carlo recipe to measure microcanon- 
ical averages. In order to test its performance, we applied it to the Ising 
ferromagnet on a 32 x 32 square lattice. The exact values of g(E) are known 
up to this lattice size, thus it is a good standard to compare our numerical 
results with. Measuring the deviations relative to the exactly known values, 
we verified a decay proportional to 1/V counts, by increasing the counter 
(counts) of averaged samples over at least 6 decades. That is why we believe 
this microcanonical simulator presents no bias besides the normal statistical 
fluctuations. For counts ~ 10 10 , we measured relative deviations near 10 -5 
for both g(E) and the specific heat peak, obtained through BHM relation. 

PACS: 75.40.Mg Numerical simulation studies 
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Monte Carlo methods are applied to statistical physics in order to measure the thermal 
average 



< Q >T - — 7 7-1 /rr s 

of some macroscopic quantity Q (magnetisation, density, etc). The temperature T is fixed, 
and the Boltzmann constant is set to unity. Both sums run over all possible microstates 
£ available for the system. The energy (quantity Q) corresponding to S is denoted by E$ 
(Qs)- Within traditional computer simulations, instead of taking all possible states, one 
takes only a finite set of them, i.e. a Markovian chain of states randomly tossed according 
to probabilities dictated by the Boltzmann exponential factors exp(— Es/T), the so-called 
importance sampling. Thus, one needs to fix a particular value for the temperature T, 
before running the computer job. In order to determine the full dependence of < Q >t 
upon T, one needs to run the job again and again, for different values of T. 
An alternative is to re-write the same average as 

„ n ^ _ Y.E < g(g) > S(E) exp(-g/T) 

<Q>T Y. E 9(E)eM-E/T) ' (2) 

where the degeneracy g(E) counts the number of states with energy E, and both sums run 
over all possible energies. The microcanonical average 

of the same quantity Q corresponds to a fixed energy, i.e. the sum in (3) runs only over the 
states S[E] belonging to energy level E. This microcanonical average (3) is simpler than 
the canonical counterpart (1) or (2), prescribing exactly the same weight to all averaging 
states, i.e. it is a uniform averaging process within each energy level separately. 

Only the Boltzmann factor exp(—E/T) appearing in equation (2) depends on the 
temperatute, carrying all thermodynamic information about the environment which con- 
tinuously exchanges energy with the system under study. Contrary to this, both g(E) and 
< Q(E) > are independent of the particular way this energy exchange occurs, independent 
of the environment: they are more fundamental properties of the system alone, defined only 
by its energy spectrum. They are not thermodynamic quantities, and do not depend upon 
thermodynamic concerns like temperature, equilibrium, etc. In practical terms, equation 
(2) allows one to determine < Q >t for any value of T, from the knowledge of the energy 
functions g (E) and < Q(E) >. 
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The Broad Histogram Method (BHM) [1] relates g(E) with the microcanonical aver- 
ages of two macroscopic quantities iV up and iV dn , measured at the current state S. First, 
one needs to adopt some protocol of allowed movements which could be performed on 
S, leading to another possible state S' . One can adopt any such a protocol, the only 
restriction being its reversibility (if S — > S' is allowed, so is S' — > S). Considering the 
Ising model, for instance, the protocol could be chosen to be the whole set of single-spin 
flips (among many other alternative choices). Given such a protocol, N up (S) counts the 
number of allowed movements which could be performed on S, leading to a fixed energy 
increment AE (which must be chosen a priori, although the method is also independent 
of this choice). Analogously, N dn (S) counts the number of allowed movements decreasing 
the energy of S by the same fixed amount AE. The fundamental BHM relation [1] is 

g(E) < N up (E) > = g(E + AE) < N dn (E + AE) > , (4) 

where the microcanonical averages of iV up and iV dn are defined by equation (3). By 
knowing these energy functions, equation (4) allows one to determine g(E) along all the 
energy axis, in steps of AE — a constant, unimportant pre-factor is cancelled out by 
performing the average in equation (2). Relation (4) is shown to be valid for any energy 
spectrum, under completely general grounds [2]. Also, the same reasoning could be applied 
for another basic quantity q, instead of the energy, by considering the degeneracies g(q) 
instead of g(E). For simplicity, we will restrict ourselves to the case of the energy. 

Thus, the Broad Histogram Method consists in: i) to choose some protocol of al- 
lowed movements, as well as an energy jump AE; ii) to measure, by any means, the 
microcanonical average < Q(E) > as a function of the energy, as well as < N up (E) > and 
< N dn (E) > which determine g(E) through equation (4); Hi) to obtain the desired thermal 
average through equation (2). There is no approximation at all, and the final numerical 
accuracy depends exclusively upon the microcanonical measuring strategy adopted in step 
ii. The method is reviewed in [3]. Some references where it is used are [4-16]. 

Let's briefly analyse hereafter some possible computer strategies one can adopt in 
order to measure the microcanonical averages < Q(E) >, < N up (E) > and < N dn (E) >, 
as functions of the energy. A Markovian chain of states is obtained by performing random 
movements transforming the current state into another. These movements are tossed 
among some previously defined set of possibilities, another protocol which has nothing to 
do with the BHM protocol of virtual movements. Both protocols could even be chosen to 
be the same, but not necessarily. 

The first direct strategy is to keep always the same fixed value E: by starting from 
some state corresponding to the desired energy, one simply rejects any tossed movement 
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which changes the current energy. After one has already a large enough number of visited 
states inside this particular energy level, the same process is repeated for other levels. De- 
pending on the adopted protocol of (real) movements, this strategy could lead to ergodicity 
problems: due to its high rejection rate, one risks to sample only a biased sub-set of states 
belonging to energy level E, violating the required uniform visitation. Moreover, it is also 
an inefficient strategy, in what concerns the computer time, again due to its high rejection 
rate. 

The opposed alternative is to accept also movements leading to energy jumps, storing 
separated averages for each energy level, in parallel. Obviously, this option is much more 
efficient in what concerns the computer time. However, one cannot simply accept any 
tossed movement: energy increments would occur more often than decrements, because 
g(E) is normally a fast increasing function of the energy. As a result, at the end, only 
states corresponding to the region near the maximum of g(E) would be sampled. Thus, 
some movement-rejection prescription must be adopted, and we get back to the uniformity 
violation problem. One can adopt some already well established movement-rejection pre- 
scription based on detailed balance arguments. For instance, canonical, fixed-temperature 
dynamics could be adopted [4,6], sampling states inside the narrow energy window corre- 
sponding to the fixed value of T. In order to get results on a broader energy window, one 
can simply superimpose the histograms obtained for different computer runs corresponding 
to different values of T. 

Another possibility is to adopt one of the many multicanonical dynamics [7,9,11,13,14]. 
Within this approach, one tunes the independent rejection rate during the computer run, 
in order to obtain the same visit probability for all energy levels, i.e. a flat distribution 
along the energy axis, at the end. In this case, the visitation probability to each particular 
state would be proportional to l/g(E). The so-called multicanonical methods [17-19] are 
based just on this feature: by recording the acceptance probability for energy-increasing 
movements E — > accumulated during the run and which must be equal to g(E) / g(E') 
at the end, one gets the function g(E) except for an unimportant global factor which 
cancels in equation (2). BHM is completely distinct from multicanonical approaches in 
many features. In particular, the infomations extracted from each Estate S belonging to 
the averaging Markovian chain are the macroscopic values of iV up (S) and N dn {S), not 
the mere one- more- visit upgrade V(E) — > V(E) + 1. That is why BHM gives more accurate 
results than multicanonical approaches, even taking into account the same Markovian chain 
of averaging states, as shown in [15] where the multicanonical dynamics [18] is adopted in 
order to measure < N up (E) > and < N dn (E) >. At the end, g{E) is determined twice, by 
following the multicanonical traditional way or, alternatively, by the BHM equation (4), 
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both using data taken from the same set of averaging states. A clear accuracy advantage for 
BHM is reported [15]. Moreover, due to the macroscopic character of the BHM quantities 
iV up and iV dn this advantage still increases for larger and larger systems. 

Another crucial difference relative to multicanonical approaches is that BHM requires 
only the uniformity of visits among the states inside each energy level, separately, in 
order to get the correct microcanonical averages. BHM does not need any detailed balance 
between visits to different levels E and E', nor the multicanonical flat distribution along 
the energy axis. Any dynamic strategy which is good for multicanonical methods will be 
also good for BHM (besides the accuracy advantage quoted in the last paragraph), but the 
reverse is not true. Thus, within BHM, other not-so-restricted dynamic strategies could 
be used. 

Profiting from this feature, we decided to test a very simple dynamic strategy inspired 
by reference [20] (although within a different method, the dynamic rule introduced in [20] 
is essentially the same as presented hereafter). The idea is to avoid movement rejections 
at all, within the energy level currently being sampled for averaging purposes. Rejections 
will be restricted to other energies, whose states are never included into the averaging 
statistics. Let's consider that the maximum energy jump allowed by the adopted protocol 
of (real) movements corresponds to n levels above or below the current energy E. Then, 
let's take an energy window of 2n + 1 adjacent levels, starting from some state inside it. 
Any movement which keeps the system still inside this window will be accepted. This is 
the dynamic strategy we propose here. Averages are taken only for the central level E 
inside the chosen energy window: the system is allowed to visit the other n levels above 
it, as well as the n levels below it, nevertheless without measuring anything during these 
side visits. Note that no tossed movement will be rejected, if the current energy is just E. 
Thus, the averaging process is completely rejection-free, avoiding any systematic bias due 
to artificial rejections rules. The same strategy can be easily applied also for continuous 
energy spectra, by taking averages only within a narrow, rejection-free energy window 
centered inside another broader, free-visit window: any movement to outside this latter 
would be rejected. 

In what follows, we consider a L x L square lattice Ising ferromagnet, with L = 32 for 
which the exact function g(E) is known [21]. The energy of the current state S is counted 
as the total number of its unsatisfied bonds, i.e. the total number of neighbouring pairs 
of spins pointing one up and the other down. The energy spectrum corresponds to all 
even numbers between and 2, 048, i.e. E = 0, 2, 4, 6, 8 . . . 2,048, which can also be 
represented by energy densities e = E/(2L 2 ). The degeneracies are g(E) = 2, 0, 2,048, 
4,096, 1,057,792 ... 2, respectively. This spectrum is symmetric in relation to its center 
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at £ max = 1,024 (or e = 0.5), where g(E max ) « 6.3 x 10 306 . We need only its first half, 
corresponding to positive temperatures. The critical energy (at the thermodynamic limit) 
corresponds to e c « 0.147 (E c m 300, for L = 32). 

We will adopt single-spin flips as the protocol of movements (both for the real move- 
ments tossed during the computer run, as well as for the virtual ones we consider in order to 
count iV up and iV dn ). Starting from the current state S with energy E, 1,024 movements 
would be available, by tossing one random spin to flip. They can be classified according to 
the possible energy jumps, E — > E ± AE, where AE could be 0, 2 or 4. Thus, our energy 
window will have 5 adjacent levels, the central one, E, where the averages will be mea- 
sured, which is rejection-free, plus two neighbouring levels above it, and two others below 
it. During the random walk performed inside this window, we measure the values of N up 
and iV dn for the current state, every time its energy is E, and accumulate them into two 
-E-histograms H up (E) and H dn (E). Also a visit counter V(E) is upgraded to V(E) + 1. If 
the energy of the current state is not E, nothing is measured or accumulated. At the end, 
we take the averages < N up (E) >= H up (E)/V(E) and < N dn (E) >= H dn (E)/V(E). 
Note that this is the only role played by the final V(E): no comparison with the neigh- 
bouring values V(E± AE) is needed, no further information is extracted from these values. 
They must only be large enough in order to provide a good statistics. For each new E, a 
new 5-levels energy window centered on it is sampled, and the whole process is repeated. 

By following this dynamic rule and by using the BHM relation (4), we measured the 
quantity ln[g(E + AE)/g(E)] for 15 adjacent levels E = 300, 302 .. . 328, at the critical 
region. The deviations from the exact values were averaged (root mean square) over these 
15 levels, and are shown as a function of the number of averaging states sampled inside 
each energy level (counts), in figure 1. The squares corresponds to AE = 4, while the 
diamonds represent AE = 2. The dashed straight line (1/V counts) indicates that no 
systematic errors besides the normal statistical fluctuations occur, giving credit to our 
simple dynamic rule. According to these results, to improve the numerical accuracy is a 
simple matter of taking more and more averaging states inside each energy level, up to the 
computer time available. 

In order to perform thermal averages, one does not need the same accuracy along the 
whole energy axis. The function [g(E) exp (— E/Tc)} 2 is displayed by the dotted line in 
figure 2, where T c = 2.293930 is the exact location of the specific heat peak. This curve 
displays the (squared) relative contribution of each energy to the partition function. As 
the number of sampled averaging states inside each energy level is proportional to the 
squared numerical accuracy, the dotted line in figure 2 shows the ideal profile of visits 
one needs in order to have equally accurate contributions from each energy level. It is a 
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sharp peaked curve, according to which the computational effort can be concentrated only 
inside a narrow energy window. Profiting from this feature, we have shaped the profile of 
visits displayed by the solid line, in figure 2. The possibility of designing this profile of 
visits, sampling different numbers of averaging states for different energies, according to 
the relative contribution of each energy region, is a further big advantage of BHM over 
multicanonical methods. Almost all the computational effort is concentrated near the 
peak. 

Figure 3 shows a detailed comparison of our simulational results with the exact spe- 
cific heat curve, near its peak. Being a derivative, which corresponds to a mathematical 
procedure in which numerical accuracy is strongly compromised, this quantity is a good 
standard for worst-case comparisons, moreover near its peak. Nevertheless, the relative 
deviations we obtained are compatible with the number (counts = 8 x 10 9 ) of averaged 
states per energy level, at the peak. Moreover, better yet accuracies were obtained along 
all the temperature range from to oo, always by using the same simulational data, i.e. 
the same averaged values for the BHM quantities iV up and iV dn measured during a single 
computer run. 

In short, we have tested a simple dynamics which is very efficient in measuring mi- 
crocanonical averages. It is essentially the same dynamics as introduced in [20], for other 
purposes. Here, the aim is to measure the microcanonical averages of some particular 
macroscopic quantities defined within the broad histogram method [1-3]. Once one knows 
these averaged quantities as functions of the energy, the method provides the energy de- 
generacy function g(E) through an exact relation. During the same computer run, the 
microcanonical average < Q(E) > of the quantity Q of interest is also measured. Then, 
once one knows g(E) and < Q(E) >, the canonical thermal average < Q >t can be 
determined by equation (2), for any, continuously varying temperature T, without 
resorting again to further computer simulations. According to our tests, the cur- 
rent dynamic rule does not introduce any systematic averaging bias, besides the normal 
statistical fluctuations which decay proportionally to 1/ y 'counts. Thus, by applying this 
dynamic rule to BHM, to improve more and more the numerical accuracy is a simple mat- 
ter of increasing the computer time. Among the further advantages of the microcanonical 
dynamics tested here, we can quote: i) its implementation simplicity, without detailed 
balance and other complications; ii) no movement-rejections at all, within the averaging 
energy level; Hi) the possibility to shape the profile of visits along the energy axis, accord- 
ing to the desired accuracy; iv) no randomness at all is used in order to decide to perform 
or not the currently tossed movement [22]; v) short and non-periodic waiting time between 
consecutive averaging states [23]. 
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Figure Captions 

1 Deviations between measured g(E) and the exact values, as a function of the number 
(counts) of sampled averaging states. The dashed line corresponds to 1/V counts. 
Data for 32 x 32 Ising ferromagnet. 

2 Profile of visits along the energy axis, which could be shaped according to the impor- 
tance of each energy contribution to the partition function (dotted line). 

3 Detail of the specific heat peak, a worst-case comparison. 
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